// SpecOffsetCorrect.cpp
// Corrects for offsets in spectral data 
// Author: Garreth Ruane


// Usage: The first command line argument is used to define the name of the input file. The output file has a '.fxt' suffix. The second command line argument is a "trusted" wavelength within the spectrum. Any number of subsequent wavelength pairs that define junctions can be passed. The regions outside of the trusted region will be scaled in order to remove offset artifacts at junctions. 
// For Example: './SpecOffsetCorrect input.txt 1050 1000 1001 1830 1831' will correct the two outside regions to match an interior region ranging from 1001 to 1830. 

#include <iostream>
#include <cmath>
#include <cstring>
#include <stdlib.h>
#include "readSpectrumFile.h"

using namespace std;

const int MAX_NUMBER_OF_JUNCTIONS = 50;

// Searchs through array, A, and finds the element nearest to the target value. 
// Returns the index of that value. 
int channelSearch( double A[], double target, int F, int L);

void fixData( double wavelengths[], double data[], double newdata[], double lowerboundi, double upperboundi, int numberOfJunctions, int numberOfChannels );

void deleteElement( double array[], int index, int N );

int main(int argc, char *argv[])
{		

	// Read Command Line Arguments
	
	// Check format
	if(argc < 5 || argc % 2 == 0)
	{
		cout << "Invalid function call" << endl;
		cout << "SpecOffsetCorrect usage: ./SpecOffsetCorrect <input file name> <trusted wavelength> <{left edge of junction} {right edge of junction}> ... " << endl;
		cout << "For Example: './SpecOffsetCorrect input.txt 1050 1000 1001 1830 1831' will correct the two outside regions to match an interior region ranging from 1001 to 1830."  << endl;
		return -1;
	}
	
	// Read in input and output filenames
	char *inputfile;	
	char outputfile[50];
	inputfile = argv[1];
	strcpy(outputfile,inputfile);
	
	// change output file extension
	char *pExt = strrchr(outputfile, '.');
	if (pExt != NULL)
		strcpy(pExt, ".fxt");
	else
		strcat(outputfile, ".fxt");
	
	cout << endl;
	cout << "Output Spectrum File: " << outputfile << endl;
	
	// The second command line argument is a wavelength in the trusted region.
	double trustedw = atof(argv[2]);
	
	int i = 3, 			// command line argument index 
	numOfJunctions = 0; // the junction counter 
	
	// arrays to hold junction points
	// The lowerpts array holds lower junction boundary.
	// The upperpts array holds upper junction boundary.
	double lowerpts[MAX_NUMBER_OF_JUNCTIONS/2], upperpts[MAX_NUMBER_OF_JUNCTIONS/2];
	
	// Read in junction points and count them too.
	while(i < argc - 1)
	{

			if( atof(argv[i]) != 0 && atof(argv[i+1]) != 0)
			{
				lowerpts[numOfJunctions] = atof(argv[i]);
				upperpts[numOfJunctions] = atof(argv[i+1]);
						
				numOfJunctions++;
				i+=2;
			}
			else
			{
				// This shouldn't happen. 
				cout << "Bad command line arguments!" << endl;
				return -1;
			}
	}

	// Count the total number of channels
	int numchans = countChannels(inputfile); 
		
	// wavelengths - w, data - d, error-e, newd - fixed data. 
	double w[numchans],d[numchans], e[numchans], newd[numchans];

	//read the file and save info.
	readFile(inputfile,w,d,e,numchans);	
	
	// Check that Input wavelengths are within bounds. 
	for( i = 0; i < numOfJunctions; i++)
	{
		if(lowerpts[i] < w[0] || lowerpts[i] > w[numchans-1])
		{
			cout << "Error: " << lowerpts[i] <<" is out of bounds." << endl;
			return -1;
		}
		if(upperpts[i] < w[0] || upperpts[i] > w[numchans-1])
		{
			cout << "Error: " << upperpts[i] <<" is out of bounds." << endl;
			return -1;
		}
	}
	
	// Initialize newd to d
	for(i = 0; i < numchans; i++)
		newd[i] = d[i];
		
	//Start Loop	
	int passnum = 1;	
	int trustedindex;	
	bool found;
	while(numOfJunctions > 0)
	{	
		cout << endl;
		cout << "Pass " << passnum << endl;
			
		trustedindex = 0;
		while( trustedw >= upperpts[trustedindex] )
			trustedindex++;
		if(trustedindex > 0)
			trustedindex--;
			
		if( trustedw <= lowerpts[0] )
		{
			fixData(w, d, newd, w[0], lowerpts[0], numOfJunctions, numchans );	
			numOfJunctions-=1;
		}
		else if( trustedw >= upperpts[numOfJunctions - 1] ) 
		{
			fixData(w, d, newd, upperpts[numOfJunctions - 1], w[numchans - 1], numOfJunctions, numchans );
			numOfJunctions-=1;
		}		
		else
		{

			fixData(w, d, newd, upperpts[trustedindex], lowerpts[trustedindex + 1], numOfJunctions, numchans );	
			
			for(i = 0; i < numchans; i++)
				d[i] = newd[i];
			deleteElement(lowerpts,trustedindex,numOfJunctions);		
			deleteElement(upperpts,trustedindex,numOfJunctions);
			numOfJunctions-=2;
		}
		deleteElement(lowerpts,trustedindex,numOfJunctions);
		deleteElement(upperpts,trustedindex,numOfJunctions);
		
		passnum++;
	}
	
	cout << "Completed in " << passnum - 1 << " pass(es)." << endl;
	cout << endl;
	
	// Print results to file. 	
	ofstream fs;
	fs.setf(ios::fixed,ios::floatfield);
	fs.precision(6);
	fs.open(outputfile);
	for(i = 0; i < numchans; i++)
		fs << w[i] << "     " << newd[i] << endl; 
	fs.close();			
	
	return 0;
}

// Searchs through array, A, and finds the element nearest to the target value. 
// Returns the index of that value. 
int channelSearch( double A[], double target, int F, int L)
{
	int i, channel = 0;
	double diff, mindiff;
	mindiff = fabs(A[F] - target);
	
	for( i = F + 1; i <= L; i++)
	{
		diff = fabs(A[i] - target);
		if(diff < mindiff)
		{
			mindiff = diff;
			channel = i;
		}
	}
	return channel;
}

// lowerbound is lower boundary of trusted region
// upperbound is upper boundary of trusted region
void fixData( double w[], double d[], double newd[], double lowerbound, double upperbound, int numberOfJunctions, int numberOfChannels )
{
		cout << "Trusted Region: " << lowerbound << " - " << upperbound << endl;
		int lowerboundi = channelSearch( w, lowerbound, 0 , numberOfChannels - 1);
		int upperboundi = channelSearch( w, upperbound, 0 , numberOfChannels - 1);
		double m1,m2,m3;
		int i;
		if(lowerboundi > 0)
		{
			m1 = (d[lowerboundi - 1] - d[lowerboundi - 2])/(w[lowerboundi - 1] - w[lowerboundi - 2]);
			m2 = (d[lowerboundi + 1] - d[lowerboundi])/(w[lowerboundi + 1] - w[lowerboundi]);
			m3 = (m1+m2)/2.0;
		
			newd[lowerboundi - 1] = m3*(w[lowerboundi - 1] - w[lowerboundi]) + d[lowerboundi];

			for(i = 0; i < lowerboundi - 1; i++)
				newd[i] = newd[lowerboundi - 1]*d[i]/d[lowerboundi - 1];
		}
		if(upperboundi < numberOfChannels - 1)
		{
			m1 = (d[upperboundi] - d[upperboundi - 1])/(w[upperboundi] - w[upperboundi - 1]);
			m2 = (d[upperboundi + 2] - d[upperboundi + 1])/(w[upperboundi + 2] - w[upperboundi + 1]);
			m3 = (m1+m2)/2.0;

			newd[upperboundi + 1] = m3*(w[upperboundi + 1] - w[upperboundi]) + d[upperboundi];
			for(i = upperboundi + 2; i < numberOfChannels; i++)
				newd[i] = newd[upperboundi + 1]*d[i]/d[upperboundi + 1];	
		}
}


void deleteElement( double array[], int index, int N )
{
	while(index < N)
	{
		array[index] = array[index + 1];
		index++;
	}
	
}

